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Abstract 



We consider the problem of the numerical approximation of the lin- 
ear controllability of waves. All our experiments are done in a bounded 
domain Q of the plane, with Dirichlet boundary conditions and internal 
control. We use a Galerkin approximation of the optimal control operator 
of the continuous model, based on the spectral theory of the Laplace op- 
erator in Q. This allows us to obtain surprisingly good illustrations of the 
main theoretical results available on the controllability of waves, and to 
formulate some questions for the future analysis of optimal control theory 
of waves. 
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1 Introduction 



This paper is devoted to the experimental study of the exact controUa- 
bihty of waves. All our experiments will be done in a bounded domain Q 
of the plane, with Dirichlet boundary conditions and with internal con- 
trol. We use the most natural approach for the numerical computation: a 
Galerkin approximation of the optimal control operator of the continuous 
model based on the spectral theory of the Laplace operator in Q. This will 
allow us to obtain surprisingly good illustrations of the main theoretical 
results available on the controllability of waves, and to formulate some 
questions for the future analysis of the optimal control theory of waves. 

The problem of controllability for linear evolution equations and sys- 
tems has a long story for which we refer to the review of D.-L. Russel in 
[Rus78] and to the book of J. L. Lions [Lio88]. Concerning controllabil- 
ity of linear waves, the main theoretical result is the so called "Geomet- 
ric Control Condition" of C. Bardos, G. Lebeau and J. Ranch [ ], 
GCC in short, which gives a (almost) necessary and sufficient condition 
for exact controllability. This is a "geometrical optics" condition on the 
behavior of optical rays inside Here, optical rays are just straight lines 
inside reflected at the boundary according to the Snell-Descartes law 
of reflection. The precise definition of optical rays near points of tangency 
with the boundary is given in the works of R. Melrose and J. Sjostrand 
[MS78] and [MS82]. For internal control, GCC asserts that waves in a 
regular bounded domain Q are exactly controllable by control functions 
supported in the closure of an open sub-domain U and acting during a 
time T, if (and only if, if one allows arbitrary small perturbations of the 
time control T and of the control domain U) 

GCC: every optical ray of length T in Q enters the sub-domain U. 

Even when GCC is satisfied, the numerical computation of the con- 
trol is not an easy task. The original approach consists first in discretizing 
the continuous model, and then in computing the control of the discrete 
system to use it as a numerical approximation of the continuous one. 
This method has been developed by R. Glowinski et al. (see [GLL90] and 
[GHL08]), and used for numerical experiments in [AL98]. However, as 
observed in the first works of R. Glowinski, interaction of waves with a 
numerical mesh produces spurious high frequency oscillations. In fact, the 
discrete model is not uniformly exactly controllable when the mesh size 
goes to zero, since the group velocity converges to zero when the solutions 
wavelength is comparable to the mesh size. In other words, the processes 
of numerical discretization and observation or control do not commute. A 
precise analysis of this lack of commutation and its impact on the com- 
putation of the control has been done by E. Zuazua in [/juau^] and [zjuauo]. 

In this paper, we shall use another approach, namely we will dis- 
cretize the optimal control of the continuous model using a projection 
of the wave equation onto the finite dimensional space spanned by the 
eigenfunctions ej of the Laplace operator in Q with Dirichlet boundary 



3 



conditions, — Ae^ = <^|ej, ej\dn=o , with cuj < cu. Here, cu will be a cutoff 
frequency at our disposal. We prove (see lemma 2 in section 2.3) that 
when GCC is satisfied, our numerical control converges when u ^ oo 
to the optimal control of the continuous model . Moreover, when GCC 
is not satisfied, we will do experiments and we will see an exponential 
blow-up in the cutoff frequency u of the norm of the discretized optimal 
control. These blow-up rates will be compared to theoretical results in 
section 2.3. 

The paper is organized as follows. 

Section 2 is devoted to the analysis of the optimal control operator for 
waves in a bounded regular domain of R^. In section 2.1, we recall the 
definition of the optimal control operator A. In section 2.2, we recall some 
known theoretical results on A: existence, regularity properties and the 
fact that it preserves the frequency localization. We also state our first 
conjecture, namely that the optimal control operator A is a microlocal 
operator. In section 2.3, we introduce the spectral Galerkin approxima- 
tion M^]^ of A, where a; is a cutoff frequency. We prove the convergence 
of M^]^ toward A when cj ^ oo, and we analyze the rate of this conver- 
gence. We also state our second conjecture on the blow-up rate of 
when GCC is not satisfied. Finally, in section 2.4, we introduce the basis 
of the energy space in which we compute the matrix of the operator Mt,uj ■ 

Section 3 is devoted to the experimental validation of our Galerkin 
approximation. In section 3.1 we introduce the 3 different domains of the 
plane for our experiments: square, disc and trapezoid. In the first two 
cases, the geodesic flow is totally integrable and the exact eigenfunctions 
and eigenvalues of the Laplace operator with Dirichlet boundary condi- 
tion are known. This is not the case for the trapezoid. In section 3.2 we 
introduce the two different choices of the control operator we use in our ex- 
periments. In the first case (non smooth case), we use x{t,x) = l[o,T]lt7. 
In the second case (smooth case), we use a suitable regularization of the 
first case (see formulas (66) and (67)). Perhaps the main contribution of 
this paper is to give an experimental evidence that the choice of a smooth 
enough control operator is the right way to get accuracy on control com- 
puting. In section 3.3, in the two cases of the square and the disc, we 
compare the exact eigenvalues with the eigenvalues computed using the 
5-points finite difference approximation of the Laplace operator. In sec- 
tion 3.4, formula (69), we define the reconstruction error of our method. 
Finally, in section 3.5, in the case of the square geometry, we compare the 
control function (we choose to reconstruct a single eigenvalue) when the 
eigenvalues and eigenvectors are computed either with finite differences or 
with exact formulas, and we study the experimental convergence of our 
numerical optimal control to the exact optimal control of the continuous 
model when the cutoff frequency goes to infinity. 

The last section, 4, presents various numerical experiments which il- 
lustrate the theoretical results of section 2, and are in support of our two 
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conjectures. In section 4.1, our experiments illuminate the fact that the 
optimal control operator preserves the frequency localization, and that 
this property is far much stronger when the control function is smooth, 
as predicted by theoretical results. In section 4.2, we use Dirac and box 
experiments to illustrate the fact that the optimal control operator shows 
a behavior very close to the behavior of a pseudo-differential operator: 
this is in support of our first conjecture. In section 4.3, we plot the re- 
construction error as a function of the cutoff frequency: this illuminates 
how the rate of convergence of the Galerkin approximation depends on 
the regularity of the control function. In section 4.4, we present various 
results on the energy of the control function. In section 4.5, we compute 
the condition number of the matrix Mt,uj for a given control domain [/, as 
a function of the control time T and the cutoff frequency uo. In particular, 
figures 31 and 32 are in support of our second conjecture on the blow-up 
rate of M^]^ when GCC is not satisfied. In section 4.6, we perform ex- 
periments in the disc when GCC is not satisfied, for two different data: 
in the first case, every optical rays of length T starting at a point where 
the data is not small enters the control domain t/, and we observe a rather 
good reconstruction error if the cutoff frequency is not too high. In the 
second case, there exists an optical ray starting at a point where the data 
is not small and which never enters the control domain, and we observe 
a very poor reconstruction at any cutoff frequency. This is a fascinating 
phenomena which has not been previously studied in theoretical works. It 
will be of major practical interest to get quantitative results on the best 
cutoff frequency which optimizes the reconstruction error (this optimal 
cutoff frequency is equal to oo when GCC is satisfied, the reconstruction 
error being equal to in that case) , and to estimate the reconstruction er- 
ror at the optimal cutoff frequency. Clearly, our experiments indicate that 
a weak Geometric Control Condition associated to the data one wants to 
reconstruct will enter in such a study. 



2 The analysis of the optimal control op- 
erator 

2.1 The optimal control operator 

Here we recall the basic facts we will need in our study of the optimal 
control operator for linear waves. For more details on the HUM method, 
we refer to the book of J.-L. Lions [Lio88]. 

In the framework of the wave equation in a bounded open subset Q of 
with boundary Dirichlet condition, and for internal control the problem 

of controllability is stated in the following way. Let T be a positive time, 

U a non void open subset of r2, and x(^5 ^) as follows: 

x(t,x)=V^(t)xo(x) (1) 

where xo is a real function on Q, such that support (xo) = U and 
Xo(a:) is continuous and positive for x G [/, '0 G C°°([0,T]) and '0(t) > 
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on ]0, T[. For a given / = (i^o, '^^i) ^ Hq(Q) x L'^(Q), the problem is to find 
a source v{t,x) G L^(0, T; L^(r2)) such that the solution of the system 

fn\u = X'^ in]0,+oo[xQ 
u\dn = 0, t > (2) 
{u\t^o,dtu\t^o) = (0,0) 

reaches the state / = (uo,ui) at time T. We first rewrite the wave oper- 
ator in (2) as a first order system. Let A be the matrix 

.A -{I ,3, 

Then A is a unbounded self-adjoint operator on if = Hq{Q) x L'^{Q), 
where the scalar product on Hq(Q) is j^VuVvdx and D{A) — {u — 
(uo.ui) G H, A{u) G H, uo\dn = 0}. Let A = ^y—AD where —Ad 
is the canonical isomorphism from Hq(Q) onto H~'^(Q). Then A is an 
isomorphism from Hq(Q) onto L^(Q). The operator B{t) given by 

is bounded on H, and one has 

The system (2) is then equivalent to 

{dt - iA)f = B{t)g, /(O) = (6) 

with / = (u^dtu), g = (A"^^;,0). For any g{t) G L^([0, oo[, if), the 
evolution equation 

{dt-iA)f = B{t)g, /(0) = (7) 

admits a unique solution / = S{g) G C°([0, +oo[, H) given by the Duhamel 
formula 

f{t)= f ^^'-'^^B{3)g(s)ds (8) 
Let T > be given. Let IZt be the reachable set at time T 

TlT = {fe H, 3g € L\[0,T],H), f = S{g){T)} (9) 

Then 71t is a linear subspace of ii, and is the set of states of the system 
that one can reach in time T, starting from rest, with the action of an 
source g filtered by the control operator B. The control problem con- 
sists in giving an accurate description of T^t, and exact controllability is 
equivalent to the equality IZt = H. Let us recall some basic facts. 

Let Ti = L^([0, T], H). Let be the closed subspace of Ti spanned by 
solutions of the adjoint evolution equation 

J'^lheH, {dt - iA*)/i = 0, h{T) = hT eH} (10) 
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Observe that, in our context, A* = A, and the function h in (10) is 
given by h{t) — e~^^^~^^^hT- Let i3* be the adjoint of the operator 
g ^ S[g)(T). Then 23* is the bounded operator from H into Ti defined 
by 

B\hT){t) = B*(t)e-^^^-*^^/iT (11) 
For any g G ^^([0, T] , if), one has, with /t = S[g){T) and h[s) = 
e~^^^~^^^hT the fundamental identity 

(/t|/it)« = r [B[s)g{s)\Ks))ds = ig\B* {hT))n (12) 

^0 

From (12), one gets easily that the following holds true 

IZt is a dense subspace of H <^=^ i3* is an injective operator (13) 

which shows that approximate controllability is equivalent to a uniqueness 
result on the adjoint equation. Moreover, one gets from (12), using the 
Riesz and closed graph theorems, that the following holds true 

IZt^H ^ 3C, < V/iGif (14) 

This is an observability inequality, and is called the observability op- 
erator. We rewrite the observability inequality (14) in a more explicit 
form 

rT 

3C, \\h\\]j<C \\B%s)e-'^^-'^'^h\\\jds \/h e H (15) 

Assuming that (15) holds true, then IZt — ii, Im(i3*) is a closed 
subspace of 7Y, and B* is an isomorphism of H onto Im(i3*). For any 
/ G ii, let C/ be the set of control functions g driving to / in time T 



C/ = |g € L2([0,T],if), / = j\^(^-'^^B\ 



s)g{s)ds } (16) 



From (12), one gets 

Cf=go + (Imi3*)^, go G Imi3* nC/ (17) 

and go = B^hr is the optimal control in the sense that 

niin{||^||^2([o,T],i/), 9 ^ Cf} is achieved at g = go (18) 

Let A : H ^ H, A(/) = hr be the control map, so that the optimal 
control go is equal to ^o(t) = S*(t)e"'^^"*^^A(/). Then A is exactly the 
inverse of the map Mt : H ^ H with 



Mt = [ m{T - t)dt f m{s)ds 
Jo Jo 



T rT 

rn(s)ds 

(19) 



m(s) = e^"^B(T-s)B*(T-s)e-^"^ 

Observe that m(s) = m*(s) is a bounded, self-adjoint, non-negative oper- 
ator on H. Exact controllability is thus equivalent to 



3C7 > 0, Mt = / ^ e^*^-*'-" f J ^2*^^ ^ ] e-^<^-*>^dt > CId (20) 
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With A = the optimal control is then given by go{t) = 5*(t)e"'^^"*^^A(/) 

and is by (5) of the form go = {X-\dtw,0) where w{t) = e-^(^-*^^A(/) 
is the solution of 

( Bw = in R X Q, w\dQ = . . 

\ (^(T,.),a,^(T,.)) = A(/) ^^'^ 

Thus, the optimal control function v in (2) is equal to ^; = X^tw^ where 
w is the solution of the dual problem (21). The operator A = M^^, with 
Mt given by (20) is called the optimal control operator. 

2.2 Theoretical results 

In this section we recall some theoretical results on the analysis of the 
optimal control operator A. We will assume here that Q is a bounded 
open subset of with smooth boundary ^Q, and that any straight line 
in R^ has only finite order of contacts with the boundary. In that case, 
optical rays are uniquely defined. See an example of such rays in figure 1. 

[Figure 1 about here.] 
Let M = X Rt. The phase space is 

^T*M = T*M \ Tqm ^ T*M U T*aM. 

The characteristic variety of the wave operator is the closed subset E 
of ^T*M of points (x, t^^,T) such that |t| = |^| when x G Q and |t| > |^| 
when X G dO.. Let ^5**^ be the set of points 

^5'*n= {(xo,eo), with l^ol = 1 if xo G n,|eo| < 1 if xo G d^} 

For pQ — (xo,^o) G and r = ±1, we shall denote by s ^ 

(7po(<s),t — ST, t),s G R the generalized bicharacteristic ray of the wave 
operator, issued from (xo,^o,t, t). For the construction of the Melrose- 
Sjostrand flow, we refer to [ "^^78]^ [M^^^"^] and to [ ], vol 3, chapter 
XXIV . Then s 1po{s) = (x(po, s), ^(po, s)) is the optical ray starting 
at xo in the direction ^o- When xq G and |^o| < 1, then the right 
(respectively left) derivative of x(po, s) at s = is equal to the unit vector 
in R"^ which projects on G T*dQ and which points inside (respectively 
outside) Q. In all other cases, x(po, s) is derivable at s = with derivative 
equal to ^o- 

We first recall the theorem of [BLR92], which gives the existence of 
the operator A: 

Theorem 1 If the geometric control condition GCC holds true, then Mt 
is an isomorphism. 

Next, we recall some new theoretical results obtained in [DL09]. For 
these results, the choice of the control function x(^5^) (1) will be es- 
sential. 

Definition 1 The control function x{t,x) = 'ip{t)xo{x) is smooth if Xo ^ 
C°^(n) and ^(t) is flat att = Oandt = T. 
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For s G M, we denote by H^(Q, A) the domain of the operator (— ADiHchiet)*' 
One has H^{Q,A) = L^{n), H^{Q,A) = ifo(n), and if (e^)^>i is an 
orthonormal basis of eigenfunctions of —A with Dirichlet boundary con- 
ditions, — Acj — uj^Cj^ < uJi < UJ2 < one has 

H%n,A) = {/ e I>'(0),/ = ^/,e„ ^c^f |/,f < ~} (22) 

The following result of [ ] says that, under the hypothesis that 
the control function x(t, x) is smooth, the optimal control operator A 
preserves the regularity: 

Theorem 2 Assume that the geometric control condition GCC holds 
true, and that the control function x{t^^) smooth. Then the optimal 
control operator A is an isomorphism of H^~^'^(Q, A) A) for all 

s > 0. 

Observe that theorem 1 is a particular case of theorem 2 with s = 0. 
In our experimental study, we will see in section 4.3 that the regularity 
of the control function x(^5 x) is not only a nice hypothesis to get theo- 
retical results. It is also very efficient to get accuracy in the numerical 
computation of the control function. In other words, the usual choice of 
the control function x(t,x) = l[o,T]lt/ is a very poor idea to compute a 
control. 

The next result says that the optimal control operator A preserves the 
frequency localization. To state this result we briefly introduce the mate- 
rial needed for the Littlewood-Paley decomposition. Let G C°°([0, oo[), 
with 0(x) = 1 for \x\ < 1/2 and = for |x| > 1. Set 'ip(x) = 

(j){x) — (j){2x). Then '0 G C(J°(R*), '0 vanishes outside [1/4, 1], and one has 

(l){s) + ^ ^(2-'^s) = 1, Vs G [0, oo[ 

k=l 

Set ^o(s) = and = il^{2~^s) for k > I. We then define the 

spectral localization operators ipk{D), /c G N, in the following way: for 
u = aj Cj , we define 

'ipk{D)u = ^ V^fc(c^j)ajej (23) 
j 

One has ipk{D) = Id and ipi(D)2pj(D) = for \i — j\ > 2. In addition, 
we introduce 

k 

Sk{D) ^Y^^J^D) = M'^'^'D), k>0 (24) 

j=0 

Obviously, the operators ipk{D) and Sk{D) acts as bounded operators on 
H = Hq X L^. The spectral localization result of [DL09] reads as follows. 
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Theorem 3 Assume that the geometric control condition GCC holds 
true, and that the control function x(^5^) smooth. There exists C > 
such that for every k ^N, the following inequality holds true 

\\MD)A-AMD)\\h <C2-'' 
\\Sk{D)A-ASk{D)\\H <C2-'' 

Theorem 3 states that the optimal control operator A, up to lower or- 
der terms, acts individually on each frequency block of the solution. For 
instance, if en is the n-th eigenvector of the orthonormal basis of L^(Q), 
if one drives the data (0,0) to (en,0) in (2) using the optimal control, 
both the solution u and control v in equation (2) will essentially live at 
frequency Un for n large. We shall do experiments on this fact in section 
4.1, and we will clearly see the impact of the regularity of the control func- 
tion x{t,x) on the accuracy of the frequency localization of the numerical 
control. 



Since by the above results the optimal control operator A preserves 
the regularity and the frequency localization, it is very natural to expect 
that A is in fact a micro-local operator, and in particular preserves the 
wave front set. For an introduction to micro- local analysis and pseudo- 
differential calculus, we refer to [TaySl] and [Hor85]. In [DL09], it is 
proved that the optimal control operator A for waves on a compact Rie- 
mannian manifold without boundary is in fact an elliptic 2x2 matrix of 
pseudo-differential operators. This is quite an easy result, since if x{t, x) 
is smooth the Egorov theorem implies that the operator Mt given by 
(20) is a 2 X 2 matrix of pseudo-differential operators. Moreover, the ge- 
ometric control condition GCC implies easily that Mt is elliptic. Since 
Mt is self adjoint, the fact that Mt is an isomorphism follows then from 
Ker{MT) = {0}, which is equivalent to the injectivity of 23*. This is 
proved in [BLR92] as a consequence of the uniqueness theorem of Calderon 
for the elliptic second order operator A. Then it follows that its inverse 
A = M^^ is an elliptic pseudo-differential operator. 

In our context, for waves in a bounded regular open subset Q of 
with boundary Dirichlet condition, the situation is far much complicated, 
since there is no Egorov theorem in the geometric setting of a manifold 
with boundary. In fact, the Melrose- Sjost rand theorem [MS78], [MS82] 
on propagation of singularities at the boundary (see also [Hor85], vol 3, 
chapter XXIV for a proof) implies that the operator Mt given by (20) 
is a microlocal operator, but this is not sufficient to imply that its in- 
verse A is micro-local. However let po = (xo,^o) ^ T*Q, |^o| = 1 be a 
point in the cotangent space such that the two optical rays defined by 
the Melrose-Sjostrand flow (see [MS78], [MS82]) s G [0,T] 7±po{s) = 
(x(±po, s), «f (±Po, s)), with zbyOo = (xo,=b^o), starting at xq in the direc- 
tions it^o have only transversal intersections with the boundary. Then it 
is not hard to show using formula (20) and the parametrix of the wave 
operator, inside Q and near transversal reflection points at the boundary 
dQ, as presented in figure 2, that Mt is microlocally at po an elliptic 2x2 
matrix of elliptic pseudo-differential operators. 
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[Figure 2 about here.] 

More precisely, let J be the isomorphism from Hi on © given 
by 

1 / A -i 



One has 2|| Jit || ^23x^2 = ll^ll?/ieL2' u(t^x) is the solution of the 

wave operator Uu — ^ with Dirichlet boundary conditions on ^Q, and 
Cauchy data at time to equal to (1x0,1^1), then one has 

An(t, .) = Acos((t — to)A)no + sin((t — to)A)iti 

SO the effect of the isomorphism J is to split the solution u{t^ x) into a sum 
of two waves with positive and negative temporal frequency. Moreover, 
one has 

Then JMtJ~^ acts as a non negative self-adjoint operator on L^, 
and is equal to 

From (29), using the parametrix of the wave operator, inside Q and near 
transversal reflection points at the boundary ^Q, and integration by parts 
to show that T is a smoothing operator, it is not difficult to get that 
JMtJ~^ is microlocally at po G T*Q a pseudo-differential operator of 
order zero with principal symbol 



(30) 



Obviously, condition (GCC) guarantees that (Jq{J Mt J~^){po) is elliptic, 
and therefore JAJ~^ will be at po a pseudo-differential operator of order 
zero with principal symbol 

.o(.A.-)(po) = 2(^^^(-'^°) (31) 

Therefore, the only difficulty in order to prove that A is a microlocal 
operator is to get a precise analysis of the structure of the operator Mr 
near rays which are tangent to the boundary. Since the set of y9 6'*^ for 
which the optical ray 7p(s) has only transversal points of intersection with 
the boundary is dense in T*Q\Tqq (see [Hor8[]), it is not surprising that 
our numerical experiments in section 4.2 (where we compute the optimal 
control associated to a Dirac mass Sxq^xq G 11), confirms the following 
conjecture: 
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Conjecture 1 Assume that the geometric control condition GCC holds 
true, that the control function x(^5^) smooth, and that the optical rays 
have no infinite order of contact with the boundary. Then A is a microlocal 
operator. 

Of course, part of the difficulty is to define correctly what is a mi- 
crolocal operator in our context. In the above conjecture, microlocal will 
implies in particular that the optimal control operator A preserves the 
wave front set. A far less precise information is to know that A is a mi- 
crolocal operator at the level of microlocal defect measures, for which we 
refer to [ ]. But this is an easy by-product of the result of N. Burq 
and G. Lebeau in [ ]. 



2.3 The spectral Galerkin method 

In this section we describe our numerical approximation of the optimal 
control operator A, and we give some theoretical results on the numerical 
approximation Mt,uj of the operator Mr given by (20), even in the case 
where the geometric control condition GCC is not satisfied. 

For any cutoff frequency a;, we denote by Uuj the orthogonal projection, 
in the Hilbert space L^(f^), on the finite dimensional linear subspace 
spanned by the eigenvectors ej for < u. By the Weyl formula, if Cd 
denotes the volume of the unit ball in R^, one has 

N{oo) = dim(L^) (27r)"'^Vol(ll)cdCj'^ (cj ^ +oo) (32) 

Obviously, Ho; acts on H = Hq x and commutes with e**"^, A and 
J. We define the Galerkin approximation Mt,uj of the operator Mt as 
the operator on x 

Mt,^ =U^MtU^ = £ e'^^-'^^U^ (^^ ^2 J) ^^U^e-'^^-'^^dt (33) 

Obviously, the matrix Mt,u; is symmetric and non negative for the Hilbert 
structure induced by H on L^j x L^,, and by (19) one has with nuj(t) — 

B*(t)n^e-^(^-*)^ 

Mt,u. = / nl{t)n^{t)dt (34) 

^0 



1 / Q+,t^ 



By (29) one has also 

JMt,ujJ~^ = 2 \ Q 

Q±,^= /o^e±-^n.x'(T-s,.)n.'eT-^ds ^^^^ 

Let us first recall two easy results. For convenience, we recall here 
the proof of these results. The first result states that the matrix Mt,uj is 
always invert ible. 

Lemma 1 For any (non zero) control function x{tjx), the matrix Mt,uj 
is invertible. 
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Proof. Let u = (m.ui) e Ll x Ll such that Mt,uj(u) = 0. By (34) 
one has 

= (iMT,u(u)\u)H = / \\nu;(t)(u)fHdt. 
Jo 

This implies nuj{t){u) — for almost all t g]0, T[. If u{t^ x) is the solution 
of the wave equation with Cauchy data (i^o,'^i) at time T, we thus get 
by (5) and (1) il;{t)xQ{x)dtu{t,x) = for t G [0,T], and since ^(t) > 
on ]0,T[ and Xo{x) > on [/, we get dtu{t,x) = on ]0,T[xU. One has 



dtu{t^x) — <^jttj sin((T — t)ix)j)ej{x) + bj cos((T — t)a;j)ej (x) 

Thus we get Eo; <a; ^j^j^ji^) — Eo; <a; ^j^j(^) = for X G ?7, which 
implies, since the eigenfunctions ej are analytic in Q, that aj = bj = for 
all j. □ 

For any cjo < cj, we define = 1 — Ho;, and we set 

||niAna;o \\h = rA(cj, CJo) /o^\ 

Since the ranges of the operators AHujq and Mt na;o are finite dimensional 
vector spaces, one has for any coq 

lima;^oo rA(cj, ujo) - . . 

lima; ^ oo rM(c^, cJo) = 

The second result states that when GCC holds true, the inverse matrix 
M^J^ converges in the proper sense to the optimal control operator A 
when the cutoff frequency oo goes to infinity. 

Lemma 2 Assume that the geometric condition GCC holds true. There 
exists c > such that the following holds true: for any given f ^ H , let 
g — A(f), fuj = Ho;/ and g^ = M^]^(fuj). Then, one has 

<c||/-/.||h + ||A(/.)-M^,L(/-)IIh (38) 

with 

hm \\AiU)-M^Uf^)\\H=0 (39) 



Proof. Since GCC holds true, there exists C > such that one has by 
(20) (Mtu\u)h > CWufn for all ue H, hence {Mt,u;u\u)h > C\\u\\ji for 
all u ^ L^xL^. Thus, with c — C~'^ , one has ||A||i/ < c and ||M^J^ ||i^ < c 
for ah UJ. Since g - g^ = A{f - f^) + A(/^) - M'^f^, (38) holds true. 
Let us prove that (39) holds true. With Ao; = Hc^AHo;, one has 

A(/..) - M^l^U^) = niA(/.) + (A. - Mri)f^ (40) 
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Set for (jjo < (jtj, fujo,uj — (Ho; — IIc^o)/- Then one has 

m^MU)\\H = ||niAn^o(/) + niA(/.o,c.)lk 

< rA{uJ,UJ0)\\f\\H + c\\fu;o,u;\\H 

On the other hand, one has 



(41) 



= M-IIu^MtU^ - n^MT)AU (42) 

Prom (42) we get 

ll(Ac.-M^,L)/c.|k<c||MT||(||A||||/.o,c.|k+rA(c.,c^o)||/^ (43) 
Thus, for all uo < we get from (41), (43), and (40) 

l|A(/c.)-M^,:,(^)||^ < (l + c||MT||)(rA(c^,c^o)||/|k + c||/.o,-lk) (44) 

and (39) follows from (37), (44) and ||/a;o,a;|k < W^^ofWn ^ when 
ujQ ^ oo. □ 

We shall now discuss two important points linked to the previous lem- 
mas. The first point is about the growth of the function 

^^WtIoWh (45) 

when cj ^ oo. This function is bounded when the geometric control 
condition GCC is satisfied. Let us recall some known results in the general 
case. For simplicity, we assume that dO. is an analytic hyper-surface of R^. 
We know from [Leb92] that, for T > Tu, where Tu — 2 sup^^Q disto (x, t/) 
is the uniqueness time, there exists A > such that 

\og\\M-HH 

lim sup_^ ^" < A (46) 

UJ 

On the other hand, when there exists po G T*Q such that the optical ray 
s G [0, T] (s) has only transversal points of intersection with the 

boundary and is such that x(po, s) ^ U for all s G [0, T], then GCC is not 
satisfied. Moreover it is proven in [Leb92], using an explicit construction 
of a wave concentrated near this optical ray, that there exists B > such 
that 

\og\\M-l\\H 

lim inf^^oo > B (47) 

uo 

Our experiments lead us to think that the following conjecture may 
be true for a "generic" choice of the control function x(t,x): 

Conjecture 2 There exists C{T,U) such that 

\og\\M-l\\H 

lim ^" = C(T, U) (48) 



UJ^OO 
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In our experiments, we have studied (see section 4.5) the behavior of 
C(T, U) as a function of T, when the geometric control condition GCC is 
satisfied for the control domain U for T > To . These experiments confirm 
the conjecture 2 when T < Tq. We have not seen any clear change in the 
behavior of the constant C (T, U) when T is smaller than the uniqueness 
time Tu. 

The second point we shall discuss is the rate of convergence of our 
Galerkin approximation. By lemma 2, and formulas (43) and (44), this 
speed of convergence is governed by the function rA(a;, a;o) defined in (36). 
The following lemma tells us that when the control function is smooth, 
the convergence in (37) is very fast. 

Lemma 3 Assume that the geometric condition GCC holds true and that 
the control function x(t, x) is smooth. Then there exists a function g with 
rapid decay such that 

rK{u,u^)<g(^^^ (49) 

Proof. By theorem 2, the operator A*AA~* is bounded on H for all 
s > 0. Thus we get, for all s > 0: 

where we have used ||A^na;o||i^ < and ||nJ;A~*||if < uj~^ . The proof 
of lemma 3 is complete. □ 

Let us recall that JMt,ujJ~^ and the operators Q±,uj and Xj are de- 
fined by formula (35). For any bounded operator M on L^, the matrix 
coefficients of M in the basis of the eigenvectors Cn are 

Mi J = {Mei\ej) (50) 

Prom (1) and (35) one has, for uii < ui, uij < ui: 

Since '0(t) G has support in [0,T], we get that, for any /c G N, there 
exists a constant Ck independent of the cutoff frequency cj, such that one 
has 

s\ip\{uji - (jOj)^Q±,uj,i,j \ < Ck (51) 

Moreover, by the results of [ ], we know that the operator T deffned 
in (29) is smoothing, and therefore we get 

sup |(cJ^ + cjj)''|71;,i,j < Cfc (52) 

Figure 3 shows the log of JMt,ujJ~^ coefficients and illustrates the 
decay estimates (51) and (52). A zoom in is done in ffgure 4, so that we 
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can observe more precisely (51). In particular we can notice that the dis- 
tribution of the coefficients along the diagonal of the matrix is not regular. 
Figure 5 presents the same zoom for JM'lJ-^ This gives an illustration 
of the matrix structure of a microlocal operator. Figure 6 represents the 
log of JMt,ujJ~^ coefficients without smoothing. And finally, figure 7 gives 
a view of the convergence of our Galerkin approximation, as it presents 
the matrix entries of J(Aa; — M^J^)J~^, illustrating lemma 2, its proof, 
and lemma 3. 

[Figure 3 about here.] 
[Figure 4 about here.] 
[Figure 5 about here.] 
[Figure 6 about here.] 
[Figure 7 about here.] 



2.4 Computation of the discrete control operator 

For any real co, let N{oo) = sup{n,cOn < oo}- Then the dimension of the 
vector space is equal to N((jj). Let us define the following {(j)j)i<j<2N{ojy- 



0^. = £l for l<i<iVH 
i,j = ej_jv(J) for N{u) + l<j< 2N{uj) 



(53) 



Then {(l)j)i<j<2N(uj) is an orthonormal basis of the Hilbert space H^j — 

In this section we compute explicitly {MT(t>i\(t>k) ^ for all 1 < /c, / < 
2N{uj). We recall 












" 









cos{suJi)ei{x) 
-LJi sin{soOi)ei{x) 

sm{suji)ei{x)/uji 
cos{scOi)ei{x) 



(54) 



We now compute the coefficients of the Mt matrix, namely Mtu,^ 

{MT(l)n\(l)m)^'- 



MTn,m = (MT0n|0m)^ 

= {e''^BB''e-^'^ct>n\4>m)^dt 



(55) 



^dt 



We now have to distinguish four cases, depending on m, n being smaller 
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or larger than N{uj). For the case {m,n) < N{uj) we have: 




= loi 
= /o"( 
= loi 



x' 






COs{sUJn)fn{x) 
Un S\n{sUJn)fn{x) 



COs{sUJm)fm{x) 
UJm Sm{sUm)fm{x) 





UJnX^ sin{sUJn)fn{x) 



COs{siOn)fm(x) 
UJm Sm{sUJm)fm{x) 



= Jo W(t)xo{x))'^UJn Sm{sUJn)fn{x)\uJm Sm{sUJm) fm{x)) ^2 ds 

= '0^ sm{sun) sm{sujm) ds Xoen{x)em{x) dx 

— 0/n,mGn,m 



where 



and 



^0 



Gn,m= / Xo{x) em{x) en{x) dx 

Similarly, for the case n > N{u;), m < N{u;) we have: 




(56) 
(57) 

(58) 



- /o^(( 

= !o{ 







x' 



X^ COs{sUn)fn{x) 



Sm{siOn)fn{x)/iOn 
COs{sC0n)fn{x) 



COs(s(jJm)fm(x) 
00m Sm{sOOm)fm{x) 



COs(sUJn)fm{x) 
Um Sm{sUJm)fm{x) 

= f^{x^ COs{sUJn)fn{x)\u;mSm{sUJm)fm{x))L2(^Q) ds 

= cos(sa;n) sin(sa;m) ds Xoen{x)em{x) dx 

bn,mGn,m 



where 

bn,m = / cos(sa;n) sin(sa;rri) ds 
Jo 

For n < N{uj) and m > N{uj) we get: 
where 

Cn,m = bm,n = I COs{sUm) sin(sCt;n) ds 

Jo 

And for m,n > N{uj): 

-^Tri,m ~ dn,mGri,m 

where 



■■ / Ip'^ COs{sLUm) COs{sLUn) ds 

Jo 



(59) 
(60) 

(61) 
(62) 

(63) 
(64) 
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The above integrals have to be implemented carefully when \ix}n — 0Jrn\ is 
small, even when '0(t) = 1. 

3 Numerical setup and validation 
3.1 Geometries and control domains 

The code we implemented allows us to choose the two-dimensional domain 
Q, as well as the control domain U . In the sequel, we will present some 
results with three different geometries: square, disc and trapezoid. For 
each geometry, we have chosen a reference shape of control domain. It 
consists of the neighborhood of two adjacent sides of the boundary (in the 
square), of a radius (in the disc), of the base side (in the trapezoid). Then 
we adjust the width of the control domain, and also its smoothness (see 
next paragraph). Figures 8, 9 and 10 present these domains, and their 
respective control domains, either non-smooth (left panels) or smooth 
(right panels). 



3.2 Time and space smoothing 

We will investigate the influence of the regularity of the function x(t, x) — 
ip{t)xo{x). Different options have been set. 

Space-smoothing. The integral (58) defining Gn,m features xo- In 
the literature we find xo = 1?7, so that 



In [DL09] the authors show that a smooth Xo leads to a more regular 
control (see also theorem 3 and lemma 3). Thus for each control domain 
U we implemented both smooth and non-smooth (constant) cases. The 
different implementations of xo are: 

• constant case: xo{x,y) = lu, 

• "smooth" case: xo{x,y) has the same support of U, the width a 
of the domain {x ^ < Xo{x) < 1} is adjustable, and on this 
domain x is a polynomial of degree 2. For example, in the square 



[Figure 8 about here.] 
[Figure 9 about here.] 
[Figure 10 about here.] 




(65) 



we have: 




(66) 
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Time-smoothing. Similarly, the time integrals (57,60,62,64) defining 
a, 6, c and d features ip{t), which is commonly chosen as 1[o,t]- As pre- 
viously, better results are expected with a smooth '0(t). In the code, the 
integrals (57,60,62,64) are computed explicitly, the different implementa- 
tions of being: 

• constant case ip = 1[o,t]5 

• "smooth case" 

^(t) = ^^^^l[o,T] (67) 

3.3 Validation of the eigenvalues computation 

The code we implemented has a wide range of geometries for Q. As it 
is a spectral-Galerkin method, it requires the accurate computation of 
eigenvalues and eigenvectors. We used Matlab eigs^ function. Figure 

11 shows the comparison between the first 200 exact eigenvalues in the 
square, and those computed by Matlab with 500 x 500 grid-points. Figure 

12 presents the same comparison in the disc, for 250 eigenvalues, the 
"exact" ones being computed as zeros of Bessel function. 

[Figure 11 about here.] 

[Figure 12 about here.] 



3.4 Reconstruction error. 

In the sequel, we will denote the input data u = (uq, ui), and its image by 
the control map w = {wo,wi) = A{uo,ui), which will often be called the 
"control". We recall from section 2.1 that for a given data u = (uq, ui) to 
be reconstructed at time T, the optimal control v{t) is given by 

v{t) = xOte-'^^-'^'^w = xdte-'^^-'^'^Aiu) (68) 

Then, solving the wave equations (2) forward, with null initial conditions 
and x'^ cLS a forcing source, we reach y = (2/0,^1) in time T. Should the 
experiment be perfect, we would have (yo^yi) = {uo,ui). The reconstruc- 
tion error is then by definition: 



^0 - j/0||^l(^) + 11^1 -^l|l'2(^) 

II^oII?/i(q) + II^iIIl2(q) 




3.5 Validation for the square geometry 

3.5.1 Finite differences versus exact eigenvalues 

In this paragraph, we compare various outputs for our spectral method, 
when the eigenvalues and eigenvectors are computed either with finite- 
differences or with exact formulas. In this first experiment, we have N x 

www.mathworks . com/access/helpdesk/help/techdoc/ref /eigs .html 
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N = 500 X 500 grid-points and we use A^e = 100 eigenvalues to compute 
the G and Mr matrices. The data {uo,ui) is as follows: 



uo - ei 
1/1 = 



(70) 



Where Cn denotes the n-th exact eigenvector. The control time T is equal 
to 3, the control domain U is 0.2 wide, and we do not use any smoothing. 
For reconstruction we use 2000 eigenvalues and eigenvectors. 
Table 1 shows the condition number of the Mt matrices, and reconstruc- 
tion errors, which are very similar for both experiments. 

[Table 1 about here.] 

Figure 13 shows the relative reconstruction error between the data u and 
the reconstructed y for both experiments: 



and similarly for ui and ^i, where [/o,n is the n-th spectral coefficient of 
the data i^o, n-ih spectral coefficient (in the basis defined 

by formulas (53) in section 2.4) of the reconstructed yo when the control 
w is obtained thanks to exact eigenvalues and is the n-th spectral 
coefficient of yo when the control w is obtained thanks to finite differences 
eigenvalues. The norm ||?7o + in our basis is given by: 



For exact eigenvalues, we can see that the errors are negligible on the 
first 100-th spectral coefficients, and quite small on the next ones. We 
have similar results for finite differences eigenvalues, except that we have 
an error on the 50-th coefficient. This error does not occur when the 
reconstruction is done with the same finite differences eigenvectors basis, 
and it can probably be explained as follows: to compute the reconstructed 
y from the finite difference control we first compute an approximation of 
It; as a function of (x, y) (i.e. on the grid) from its spectral coefficients (on 
the finite differences eigenvectors basis) , then we compute the coefficients 
of this function on the exact basis (thanks to a very simple integration 
formula). We thus introduce two sources of errors, projection on the grid 
and projection on the exact basis, which do not have anything to do with 
our spectral Galerkin method. Therefore we will not discuss the matter 
in further detail here. 



Relative reconstruction error FD^ 



Relative reconstruction error SPn 



\\Uo + Ui\\ 
\\Uo + Ui\\ 



(71) 




(72) 



[Figure 13 about here.] 
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3.5.2 Impact of the number of eigenvalues 

In this paragraph, we still use the same data (70), but the number of eigen- 
values and eigenvectors A^e used to compute the Mr matrices is varying. 
Table 2 shows Mt condition numbers and reconstruction errors for various 
A^e with exact or finite-differences-computed eigenvalues. The reconstruc- 
tion is still performed with 2000 exact eigenvalues. We can see that the 
finite differences eigenvalues lead to almost as good results as exact eigen- 
values. We also observe in both cases the decrease of the reconstruction 
error with an increasing number of eigenvalues, as predicted in lemma 3. 
A 5% error is obtained with 70 eigenvalues (the input data being the 50-th 
eigenvalue), and 100 eigenvalues lead to less than 2%. 

[Table 2 about here.] 

4 Numerical experiments 
4.1 Frequency localization 

In this subsection, the geometry (square) as well as the number of eigen- 
values used (200 for HUM, 2000 for verification) are fixed. Note also that 
in this paragraph we use only exact eigenvalues for HUM and verification. 
The data is also fixed to a given eigenmode, that is: 

Uo ^ 650 Ui =0 (73) 

where en is the n-th eigenvector of —A on the square. 

[Figure 14 about here.] 

The first output of interest is the spreading of w spectral coefficients, 
compared to u. Figure 15 shows the spectral coefficients of the input 
(uo^ui) and the control {wo,wi) with and without smoothing. As pre- 
dicted by theorem 3 and lemma 3 we can see that the main coefficient 
of (wo,wi) is the 50-th of wq, and also that the smoothing noticeably 
improves the localization of w. 

[Figure 15 about here.] 

Similarly we can look at the spectral coefficients of the reconstruction 
error. Figure 16 presents the reconstruction error (see paragraph 3.4 for 
a definition) with or without smoothing. We notice that the errors occur 
mostly above the cutoff frequency (used for Mt,uj computation, and thus 
for the control computation) . Another important remark should be made 
here: the smoothing has a spectacular impact on the frequency localiza- 
tion of the error, as well as on the absolute value of the error (maximum 
of 2.10"^ without smoothing, and 8.10"^ with smoothing), as announced 
in theorem 3 and lemma 3. 

[Figure 16 about here.] 
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Remark 1 For other domains, such as the disc and trapezoid, as well 
as other one-mode input data, we obtain similar results. The results also 
remain the same if we permute uO and ul, i.e. if we choose i^O = and 
ul equal to one fixed mode. 

4.2 Space localization 
4.2.1 Dirac experiments 

In this section we investigate the locaUzation in space. To do so, we use 
"Dirac" functions ^(x,y)=(xo,?/o) data, or more precisely truncations to 
a given cutoff frequency of Dirac functions: 



where Ni is the index corresponding to the chosen cutoff frequency, with 
Ni — 100 or 120 in the sequel. Figure 17 shows the data uq and the control 
wq in the square with exact eigenvalues, without smoothing, the results 
being similar with smoothing. We can see that the support of wq is very 
similar to uqs. Figure 18 presents the reconstruction error associated to 
this experiment. We can see as before that the smoothing produces highly 
reduced errors. 



Similarly, we performed experiments with numerical approximation of 
a Dirac function as input data in the disc and in a trapezoid. Figures 
19 and 20 present the space-localization of uq and wq without smoothing 
(we get similar results with smoothing). As previously, the control wq is 
supported by roughly the same area than the input uq. In the disc we 
can see a small disturbance, located in the symmetric area of the support 
of uq with respect to the control domain U . However, this error does not 
increase with Ni^ as we can see in figure 21 (case Ni — 200) so it remains 
compatible with conjecture 1. 

Figure 22 shows the reconstruction errors for these experiments, with or 
without smoothing. As before we notice the high improvement produced 
by the smoothing. We get similar errors in the trapezoid. 



4.2.2 Box experiments in the square 

In this paragraph we consider the case uq = Ibox, where box — [0.6, 0.8] x 
[0.2,0.4] is a box in the square. The control domain U is 0.1 wide: U — 
{x < 0.1 and y > 0.9}. These experiments were performed in the square 



^0 = Z^^Ji en(xo,2/o) en 
Ul = 



(74) 



[Figure 17 about here.] 
[Figure 18 about here.] 



[Figure 19 about here.] 
[Figure 20 about here.] 
[Figure 21 about here.] 
[Figure 22 about here.] 
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with 1000 exact eigenvalues for the Mt matrix computation, the input 
data uo being defined thanks to 800 eigenvalues. Figures 23 and 24 show 
the space localization of the data uq and the control wq without and 
with smoothing. As before we can notice that the space localization is 
preserved, and that with smoothing the support of wq is more sharply 
defined. Figures 25 and 26 show the reconstruction errors for two different 
data, the first being the same as in figure 23, and the second being similar 
but rotated by 7r/4. We show here only the case with smoothing, the 
errors being larger but similarly shaped without. We can notice that the 
errors lows and highs are located on a lattice whose axes are parallel to 
the box sides. This is compatible with the structure of the wave-front set 
associated to both input data. 

[Figure 23 about here.] 

[Figure 24 about here.] 

[Figure 25 about here.] 

[Figure 26 about here.] 

4.3 Reconstruction error 

In this section we investigate lemma 3 or more precisely the subsequent 
remark ??. This remark states that the reconstruction error should de- 
crease as the inverse of the cutoff frequency without smoothing, and as 
the inverse of the fifth power of the cutoff frequency with smoothing. To 
investigate this, we perform a "one- mode" experiment (see paragraph 4.1) 
using the 50-th mode as input data. We then compute the control with an 
increasing cutoff frequency, up to 47 (finite differences case) or 82 (exact 
case), and we compute the reconstruction error, thanks to a larger cutoff 
frequency (52 in the finite differences case, or 160 in the exact case). 
Figure 27 represents the reconstruction error (with exact or finite differ- 
ences eigenvalues) as a function of the cutoff frequency (i.e., the largest 
eigenvalue used for the control function computation) . Figure 28 presents 
the same results (with finite differences eigenvalues only) for two differ- 
ent geometries: the square, and the trapezoid (general domain). The log 
scale allows us to see that the error actually decreases as the inverse of the 
cutoff frequency without smoothing, and as the inverse of the fifth power 
of the cutoff frequency with smoothing, according to remark ??. 

[Figure 27 about here.] 
[Figure 28 about here.] 

4.4 Energy of the control function 

In this paragraph we investigate the impact of the smoothing, the width 
of the control domain U and the control time T on various outputs such 
as the condition number of Mt, the reconstruction error — and the 
norm of the control function \\w\\. 

To do so we performed several one- mode experiments (see paragraph 4.1, 
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mode 500) in the square, with exact eigenvalues, 1000 eigenvalues used 
for computation of Mt, 2000 eigenvalues used for reconstruction and ver- 
ification. We chose various times: 2.5 and 8, plus their "smoothed" coun- 
terparts, according to the empirical formula Tsmooth = 15/8 * T. This 
increase of Tgmooth is justified on the theoretical level by formulas (31) 
and (30) which show that the efficiency of the control is related to a mean 
value of x(^7^) the trajectories. Similarly, we chose various width of 
U: 1/10 and 3/10, plus their "smoothed" counterpart, which are double. 
Table 3 presents the numerical results for these experiments. This table 
draws several remarks. First, the condition number of Mt, the recon- 
struction error and the norm of the control w decrease with increasing 
time and U. Second, if we compare each non-smooth experiment with its 
"smoothed" counterpart (the comparison is of course approximate, since 
the "smoothed" time and width formulas are only reasonable approxima- 
tions), the condition number seems similar, as well as the norm of the 
control function w, whereas the reconstruction error is far smaller with 
smoothing than without. 



[Table 3 about here.] 

Figures 29 and 30 emphasize the impact of the control time, they present 
the reconstruction error, the norm of the control, and the condition num- 
ber of Mt, as a function of the control time (varying between 2.5 and 
16), with or without smoothing. Conclusions are similar to the table 
conclusions. 

[Figure 29 about here.] 
[Figure 30 about here.] 



4.5 Condition number 

In this section, we investigate conjecture 2. To do so, we compute the 
condition number of Mt,uj, as we have: 

cond(MT,a;) = 1 1 Mt,c. 1 1 • 1 1 M^^^, 1 1 ^ 1 1 Mt 1 1 • 1 1 M^^^, 1 1 (75) 

Figure 31 shows the condition number of the Mt matrix as a function 
of the control time or of the last eigenvalue used for the control function 
computation. According to conjecture 2, we obtain lines of the type 

log (cond(MT,a;)) = u;.C{T, U) (76) 

Figure 32 shows for various eigenvalues numbers the following curves: 

rj. ^ log (cond(MT,a;)) ^^^^ 

Similarly, we can draw conclusions compatible with conjecture 2, as these 
curves seems to converge when the number of eigenvalues grows to infinity. 

[Figure 31 about here.] 

[Figure 32 about here.] 
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4.6 Non-controlling domains 

In this section we investigate two special experiments with non-controHing 
domains, i.e. such that the geometric control condition is not satisfied 
whatever the control time. 

First we consider the domain presented in Figure 33. 

[Figure 33 about here.] 

For this domain the condition number of the Mt matrix is large, and 
subsequently we should be experiencing difficulties to reconstruct the data 
u. We perform one-mode experiments with two different eigenvectors, one 
being localized in the center of the disc (eigenvalue 60), the other being 
localized around the boundary (eigenvalue 53) as can be seen on Figure 
34. 

[Figure 34 about here.] 

The various outputs are presented in Table 4, and we can see that the 
inversion is fairly accurate for the 53rd eigenmode, while it is logically 
poor for the 60th eigenmode. Moreover, the energy needed for the control 
process, i.e. the norm of the control is small for the 53rd eigenvector, 
while it is large for the 60-th. We can also notice that the smoothing has 
the noticeable effect to decrease the reconstruction error, the norm of the 
control function w being similar. 

[Table 4 about here.] 

In the second experiment we change the point of view: instead of 
considering one given domain and two different data, we consider one 
given data, and two different non-controlling domains. The data is again 
1^53 (see Figure 34), which is localized at the boundary of the disc. The 
first domain is the previous one (see Figure 33), the second domain is 
presented in Figure 35, it is localized at the center of the disc. 

[Figure 35 about here.] 

In either case, the condition number of the Mt matrix is large, and the 
data should prove difficult to reconstruct. Table 5 present the outputs 
we get for the two domains. As previously, we observe that the control 
process works fairly well for the appropriate control domain, with a small 
error as well as a small energy for the control. Conversely, when the 
control domain does not "see" the input data, the results are poorer: the 
energy needed is large with or without smoothing, the error is also large 
without smoothing, it is however small with smoothing. 

[Table 5 about here.] 
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Figure 1: Example of optical rays. 



32 



Figure 2: Example of optical ray with only transversal reflection points. 
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log J*M*inv(J) 
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Figure 4: View of the logarithm of the coefficients of the matrix JM^J"^, 
for the square geometry, with smooth control (zoom). This illustrates decay 
estimate (51). 
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Figure 5: View of the logarithm of the coefficients of the matrix JMj.^ J ^, for 
the square geometry, with smooth control (zoom). 
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Figure 6: View of the logarithm of the coefficients of the matrix JMtJ~^ for 
the square geometry, with non-smooth control. Note that the color scaling is 
the same as in Figure 3. 
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log abs (trunc(inv JMJ) - inv(trunc JMJ)) 
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Figure 7: View of the logarithm of the coefficients of the matrix 
J [((Mt)-1)c. - ((Mt)^)-i] = J [A^ - ((Mt)..)-!] J-\ for the square ge- 
ometry, with smooth control. The Mt matrix is computed with 2000 eigenval- 
ues, the cutoff frequency uo being associated with the 500th eigenvalue. 
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40 




Figure 10: Domain and example of a control domain for the trapezoid, with 
smoothing in space (right panel) or without (left panel). 
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Figure 11: Verification of the eigenvalues computation in the square: exact and 
finite-diff"erences-computed eigenvalues (left panel) , and their absolute diff"erence 
(right panel). 
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Figure 12: Verification of the eigenvalues computation in the disc: eigenvalues 
computed either as zeros of Bessel functions or with finite differences (left panel) , 
and their absolute difference (right panel). 



43 



0.01 2 p 
0.01 - 
_ 0.008- 
! 0.006 - 
■ 0.004- 

0.002 - 



Rel. error UO - YO SP 




0.016 
0.014 
0.012 
S 0.01 
I 0.008 
£ 0.006 
0.004 
0.002 



Rel. error UO - YO FD 




10 20 30 40 50 60 
frequency 



0.014p 
0.012 
0.01 - 

; 0.008 

i 0.006- 
0.004 - 

0.002 - 



Rel. error U1 - Y1 SP 



Rel. error U1 - Y1 FD 





Figure 13: Validation experiments in the square: Relative errors between the 
spectral coefficients of the original function u and the reconstructed function 
computed with exact eigenvalues (left panels) or finite differences eigenvalues 
(right panels), for uq and yo (top panels) or ui and yi (bottom panels). The 
errors are plotted as a function of the frequency of the eigenvalues. The com- 
putation are performed with 100 eigenvalues, corresponding to a frequency of 
about 38, the reconstruction with 2000, corresponding to a frequency of about 
160. For the readability of the figure, we plot only the major counterparts of 
the error, i.e. we stop the plot after frequency 63 (300th eigenvalue). 
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Eigenvector on the gric 




Figure 14: Representation on the grid in 3D (left panel) or contour plot (rig 
panel) of the 50-th eigenvector in the square. 
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Figure 15: One- mode experiment in the square: localization of the Fourier 
frequences of {uq^ui) (dashed line) and {wq^wi) (solid line) for a given time 
T and a given domain U without smoothing (top) and with time- and space- 
smoothing (bottom). The x-coordinate represents the eigenvalues. The input 
data uq is equal to the 50-th eigenvector, equal to an eigenvalue of about 26.8, 
and ui = 0. 
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Figure 16: One- mode experiment in the square: localization of the Fourier 
coefficients of {uq — yo^ui — ^i), where u is the data and y is the reconstructed 
function obtained from the control function tt;, for a given time T and a given 
domain U without smoothing (top panels) and with time- and space-smoothing 
(bottom panels). 
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Figure 17: Space localization of the data uq (top panels) and the control wq 
(bottom panels), for a Dirac experiment in the square, with exact eigenvalues. 
These plots correspond to an experiment without smoothing, but it is similar 
with smoothing. Left panels represent 3D view, and right panels show contour 
plots. 
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Figure 18: Difference between the data uq and the reconstructed function yo 
without smoothing (top panels) and with smoothing (bottom panels) for a dirac 
experiment in the square, with exact eigenvalues. Left panels represent 3D view, 
and right panels show contour plots. 
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Data uO Contour uO 




Figure 19: Space localization of the data uq (top panels) and the control wq 
(bottom panels), for a dirac experiment in the disc. These plots correspond to 
an experiment without smoothing, but it is similar with smoothing. Left panels 
represent 3D view, and right panels show contour plots. In this experiment, the 
input data is defined with A^^ = 100 eigenvectors. 
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Figure 20: Space localization of the data uq (top panels) and the control wq 
(bottom panels), for a dirac experiment in the trapezoid. These plots corre- 
spond to an experiment without smoothing, but it is similar with smoothing. 
Left panels represent 3D view, and right panels show contour plots. In this 
experiment, the input data is defined with Ni = 120 eigenvectors. 
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Figure 21: Space localization of the data uq (top panels) and the control wq 
(bottom panels), for a dirac experiment in the disc. These plots correspond to 
an experiment with smoothing, and it is similar without smoothing. Left panels 
represent 3D view, and right panels show contour plots. In this experiment, the 
input data is defined with Ni = 200 eigenvectors. 
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Figure 22: Difference between the data uq and the reconstructed function yo for 
a dirac experiment in the disc without smoothing (top panels) and with time- 
and space-smoothing (bottom panels). Left panels represent 3D view, and right 
panels show contour plots. In this experiment, the input data is defined with 
Ni = 100 eigenvectors. 
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Figure 23: Space localization of the control function wq (bottom panels) with 
respect to the data uq (top panels), in the square, without smoothing: 3D plots 
on the left, and contour plots on the right. 
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Figure 24: Space localization of the control function wq (bottom panels) with 
respect to the data uq (top panels), in the square, with smoothing. Left panels 
represent 3D view, and right panels show contour plots. 
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Figure 25: Difference between the data uq and the reconstructed function yo 
(top panels) and ui and yi (bottom panels) with smoothing in the square. The 
data is the identity function of a square whose edges are parallel to the x and y 
axes. Left panels represent 3D view, and right panels show contour plots. 
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Figure 26: Difference between the data uq and the reconstructed function yo 
(top panels) and ui and yi (bottom panels) with smoothing in the square. 
The data is the identity function of a square whose edges are parallel to the 
diagonals of the square. Left panels represent 3D view, and right panels show 
contour plots. 
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Figure 27: Reconstruction errors for the finite differences and exact methods, 
as a function of the cutoff frequency (i.e., the largest eigenvalue used for the 
control computation), with or without time- and space-smoothing. 
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FINITE DIFF: Reconstruction error as a function of the cutoff frequency 
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Figure 28: Reconstruction errors for the finite differences method, in the square 
and in the trapezoid, as a function of the cutoff frequency (i.e., the largest 
eigenvalue used for the control computation), with or without time- and space- 
smoothing. 
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Figure 29: Experiments in the square, with exact eigenvalues: impact of the 
smoothing on the reconstruction error (top) and on the norm of the control 
function (bottom), as a function of the control time. 
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Figure 30: Experiments in the square, with exact eigenvalues: impact of the 
smoothing on the condition number of the Mt matrix, as a function of the 
control time. 
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Condition number of HUM matrix as a function of the cutoff frequency 
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Figure 31: Condition number of the Mt^uj matrix as a function of the cutoff 
frequency uj for various control times. 
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Figure 32: Ratio of tlie log of the condition number of the Mt^cj matrix and 
the cutoff frequency cj, as a function of the control time for various eigenvalues 
numbers. 
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Figure 33: Non-controlling domain U without (left) or with (right) smoothing. 
This domain consists of the neighborhood of a radius which is truncated around 
the disc boundary. 
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Figure 34: Special modes chosen for experiment with non-controhing 
corresponding to the 53rd and 60th eigenvalues. 
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Eigenvalues computation 


Condition Number 


Reconstruction error 


Finite differences 


7.4 


1.8% 


Exact 


7.5 


1.6% 



Table 1: Validation experiments in the square: Condition numbers and valida- 
tion errors for a 100-eigenvalues-experiment in the square (without smoothing), 
where the eigenvalues are exact or computed thanks to finite differences. 
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Condition Number 
Exact Finite differences 


Reconstruction error 
Exact Finite differences 


52 


6.5 


6.4 


29.2% 


29.0% 


55 


6.6 


6.6 


17.7% 


17.6% 


60 


6.6 


6.6 


17.3% 


17.2% 


70 


7.1 


7.0 


4.9% 


5.0% 


80 


7.3 


7.2 


3.2 % 


3.3 % 


100 


7.5 


7.4 


1.6 % 


1.8 % 


200 


8.3 


8.3 


0.5 % 


1.0% 


300 


8.8 


0.3 % 


500 


9.5 


0.2 % 



Table 2: Validation experiments in the square: Condition numbers and valida- 
tion errors for various numbers A^e of eigenvalues used to compute the control 
function, where the eigenvalues are exact or computed thanks to finite differ- 
ences. The input data is the 50-th eigenvalue, and the reconstruction is per- 
formed with 2000 eigenvalues. 
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Smooth 


Width of U 


Time 


Condition number 


Rec. error 


ll^ll 


no 


1/10 


2.5 


48.8303 


0.00518843 


504.287 


yes 


2/10 


2.5 


204.048 


0.00140275 


1837.12 


yes 


2/10 


4.7 


21.7869 


4.65892E-07 


364.013 


no 


1/10 


8 


21.1003 


0.00162583 


120.744 


yes 


2/10 


8 


16.5497 


8.51442E-08 


189.361 


yes 


2/10 


15 


12.017 


8.39923E-09 


100.616 


no 


3/10 


2.5 


4.20741 


0.0014823 


147.009 


yes 


6/10 


2.5 


9.05136 


1.94519E-06 


336.704 


yes 


6/10 


4.7 


3.09927 


2.99855E-08 


125.481 


no 


3/10 


8 


3.20921 


0.000488423 


39.9988 


yes 


6/10 


8 


2.74172 


6.0204E-09 


69.8206 


yes 


6/10 


15 


2.4113 


8.55119E-10 


37.1463 



Table 3: Impact of the control time, the width of U and the smoothing, on the 
condition number of M^, on the reconstruction error and on the norm of the 
control function. These results come from one-mode experiments in the square, 
with exact eigenvalues. 
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smooth? 


Condition 
number 


53rd Eigenmode 
Rec. error \\w\\ 


60th Eigenmode 
Rec. error \\w\\ 


no 


4960 


0.63% 15 


26% 4583 


yes 


6042 


6.5 10-"^% 24 


0.11% 8872 



Table 4: Influence of the shape of the data on the reconstruction error and 
the norm of the control, with a non- controlling domain U. On the left, the 
53rd eigenmode is localized around the boundary of the circle, as is the control 
domain. On the right, the 60th eigenmode is localized around the center of the 
circle. 
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smooth? 


First domain 
Cond. nb. Rec. error \\w\\ 


Second domain 
Cond. nb. Rec. error \\w\\ 


no 


4960 0.63% 15 


3.6 10^ 68% 1.9 10"^ 


yes 


6042 6.5 10-^% 24 


3.3 10^ 9.4 10-^% 6.5 10^ 



Table 5: Reconstruction error and the norm of the control with a data 1^53 
localized at the boundary of the disc and two different non-controlling domains 
C/, the first one being localized around the boundary, the second one around the 
center. 
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